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ABSTRACT 


/ 

"-  Discrete  event  simulation*  art  computer  models  of  complex 
systems.  The  accuracy  of  the  model  In  reflecting  the  behavior 
of  the  real  system  is  determined  by,  among  other  things,  the 
values  of  parameters  for  distributions  used  in  the  model.  The 
modeller  could  allocate  experimental  resources  more  effectively 
without  loss  of  accuracy  in  the  model  if  he  or  she  could 
identify  those  parameters  to  which  the  response  of  interest 
has  the  greatest  sensitivity. 

One  method  of  doing  this  is  to  try  to  model  the  response 
as  a  polynomial  function  of  the  model  parameters.  We  are  then 
interested  in  those  terms  in  the  polynomial  which  have  non¬ 
zero  coefficients.- 

Schruben  and  Cogliano  developed  a  method  whereby  such 
polynomial  models  can  be  identified  in  fewer  computer  runs 
than  previous  methods  allowed.  Their  concept  was  to  do 
analysis  in  the  frequency  domain  rather  than  the  time  domain. 
The  virtual  independence  of  frequency  estimators  in  a  spectrum 
means  that  many  parameters  can  be  tested  independently  within 
a  single  experiment  using  spectral  methods. 

;  This  report  attempts  to  extend  the  Schruben/Cogliano 
methodology  to  cover  a  more  general  class  of  models  which 
includes  discrete-valued  parameters,  such  as  policy  decisions 
or  capacities  of  queues.  We  evaluated  the  use  of  discrete¬ 
valued  functions  as  a  basis  for  spectral  analysis.  Several 
function  sets  were  considered  as  possibilities,  and  Walsh 
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function*  were  selected  a*  th*  b**t  cholca. 

Our  preliainary  rasulta  indicata  that  Wal*h  analyai*  may 
prasant  a  proaieing  aethod  for  identifying  aignificant 
parameters.  However,  the  aethod  exhibit*  undesirable  behavior 
when  tiae  lags  are  present  in  the  aodel. 
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Section  1 


The  identification  of  parameters  or  factors  which  affect 
performance  is  an  important  area  of  operations  research  and 
statistics.  Until  recently  it  has  been  too  expensive  for  many 
computer  simulation  studies,  which  Hillier  and  Lieberman  state 
is  one  of  the  major  shortcomings  of  simulation^.  The 
traditional  approach  to  this  problem  is  to  make  separate 
computer  runs  for  each  of  many  different  factor  values.  The 
basic  experimental  unit  is  the  computer  run.  Although  the 
number  of  runs  required  may  be  reduced  using  screening 
designs^10},  it  increases  geometrically  with  the  number  of 
factors.  This  approach  becomes  prohibitively  expensive  in 
terms  of  both  user  and  computer  time  for  all  but  the  simplest 
of  models. 

Schruben  and  Cogliano^12^  (S/C)  addressed  this  problem 
utilizing  an  entirely  different  approach.  Normally  one  views 
a  parameter  as  a  fixed,  possibly  unknown,  attribute  of  the 
system.  However,  in  a  computer  simulation  the  experimenter 
has  complete  control  of  the  model  and  can  alter  parameter 
values  during  the  run.  Hence  the  terms  parameter  and  factor 
can  be  used  interchangeably  in  simulation.  Schruben  and 
Cogliano  proposed  varying  the  parameters  sinusoidally  during 
run.  Each  parameter  is  assigned  a  unique  frequency,  and 
spectral  estimators  are  used  to  analyze  the  system  output  at 
the  different  frequencies.  After  performing  a  suitable 
statistical  test,  if  the  power  spectrum  is  not  significantly 


1 


2 


different  from  tero  at  a  glvan  frequency  it  la  concluded  that 
tha  system  la  inaensitlva  to  tha  parameter  which  waa  aaaigned 
that  frequency.  By  analyzing  the  spectrum  instead  of  just  the 
assigned  frequencies,  one  can  detect  non-linear  response 
through  a  relatively  simple  set  of  relations  specified  in  the 
S/C  paper.  Most  of  their  paper  addresses  the  technical 
aspects  of  implementing  such  a  procedure. 

The  advantage  of  the  S/C  approach  is  that  analysis  is 
moved  to  the  frequency  domain.  The  output  time  series  is 
represented  using  trigonometric  functions  as  a  linear 
algebraic  basis.  The  experimental  unit  becomes  a  frequency 
band,  and  a  single  run  of  the  simulation  contains  many  almost 
independent  frequency  bands.  The  number  of  runs  required  is 
greatly  reduced. 

One  limitation  of  the  S/C  procedure  is  that  it  can  only 
be  used  to  evaluate  continuous  parameters.  This  problem  can 
be  removed  by  choosing  a  different  set  of  functions  as  a 
basis.  Three  alternative  discrete-valued  function  sets  are 
Rademacher,  Walsh,  and  Haar  functions,  which  are  named  after 
their  inventors.  It  is  the  goal  of  this  report  to  extend  the 
S/C  procedure  by  considering  the  use  of  these  alternative 
bases  for  representing  the  time  series. 

Although  this  report  concentrates  on  the  application  of 
Walsh  functions  to  computer  simulations,  the  methodology 
outlined  here  is  equally  applicable  for  any  type  of  2n 
factorial  experimental  design  situation. 


The  following  ••ctlon  gives  a  summary  of  pertlnant 
information  availabla  in  tha  litaratura  about  diacrata 
functions.  The  intention  ia  to  provide  the  reader  with  a 
aufficient  background  to  understand  the  material  in  subsequent 
sections.  The  method  proposed  in  section  2.6  for  detecting 
system  gain  apparently  has  not  appeared  in  previous 
literature.  Section  3  contains  new  material  about  the 
statistical  properties  of  the  Walsh  spectrum  estimator,  and 
tells  how  to  design  an  experiment  to  identify  significant 
parameters  in  a  model.  Section  4  gives  an  in  depth 
description  and  an  example  of  how  to  use  the  proposed 
methodology.  In  section  5  we  propose  a  number  of  extensions 
to  this  work. 


Section  2 


Our  notation  for  trigonometric  functions  will  be 

consistent  with  that  used  in  the  recent  Ph.D.  report  by  V.  J. 
f  5  1 

Cogliano1  J .  The  symbol  u  will  denote  cycles/time. 

We  will  use  the  following  definition  of  orthogonality 
(Beauchamp t 2 ^ ) . 

Defn  -  The  series  (Sn(t)}  and  (Sffi(t) )  are  orthogonal  on  the 
interval  [0,T]  if 

?  J%(t> «■<»>  *  -  { o  5 ;  i  ».» e  z+ 

and  are  orthonormal  if  k  ■  l.  (Z+  is  the  set  of  nonnegative 
integers . ) 

Any  time  series  y(t)  can  be  approximated  over  a  finite 
interval  [0,T]  by  a  weighted  sum  of  terms  in  an  orthogonal 


series: 


N-l 

y(t)  *  l  cnsn(t) 
n=0  n  n 

Sn(t)  is  term  n  of  the  orthogonal  series 
cn  is  a  weight. 

The  most  common  measure  of  precision  in  such  an 
approximation  is  the  mean  squared  error,  abbreviated  MSE. 


-T  r  N-l  -.2 

MSE  =*  y(t)  -  l  c  S  (t)  dt 
J0  L  n-0  n  n  J 


MSE  can  be  minimized  by  setting 


£ 
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T 

cn-  f  I/'1)  V*)  dt' 

It  can  be  seen  that  if  ei2nut  i,  substituted  for  6  (t)  then 

n 

this  is  the  familiar  formula  for  calculating  Fourier 
coefficients. 

A  basis  series  should  have  the  property  of 
completeness^.  One  definition  of  completeness  is  that  there 
should  not  exist  any  function  which  is  orthogonal  to  every 

element  of  the  basis.  This  implies 

lim  HSE  «  0 
N-»® 

for  all  y(t)  such  that  y(t)  contains  at  most  a  countable 
number  of  discontinuities. 

The  following  sections  present  a  brief  summary  of  the 
properties  of  three  sets  of  discrete-valued  functions  which 
might  be  considered  as  algebraic  bases  for  y(t) .  Rademacher, 
Walsh,  and  Haar  functions  are  discussed  extensively  in  all  the 
literature  for  discrete  spectral  methods.  It  turns  out  that 
only  Walsh  functions  have  all  the  properties  needed  for  the 
type  of  analysis  of  interest  to  us. 

2 . 2  Rademacher  functions 

For  the  purpose  of  evaluating  a  simulation  which  contains 
discrete-valued  parameters  we  want  to  utilize  discrete-valued 
functions  as  our  input  series.  One  such  series  is  the  set  of 
Rademacher  functions,  which  are  illustrated  in  figure  1.  These 
are  block  pulses  which  alternate  regularly  between  1  and  -l. 
Each  function  R(k,t)  is  a  function  of  the  continuous  index 
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variable  t  c  [0,T].  By  convention  T  ia  assumed  to  be  1  unless 
otherwise  stated,  but  any  value  can  be  used  with  appropriate 
scaling.  The  first  Rademacher  function,  R(0,t)  (0  £  t  &  T) , 
is  just  a  line  with  value  1.  Every  subsequent  function  is 
constructed  by  changing  the  sign  on  half  of  each  interval  in 
the  previous  function,  so  that  the  function  changes  sign  in  a 
regular  manner. 


R(0,t) 

R(l,t) 

R(2,t) 

R(3,t) 


l _ 

- I _ I - 1 _ 

1 _ I - 1 _ I - 1 _ I - 1 


+  1 
-1 

+  1 
-1 

+1 

-1 

-1 

T - 

T 


Figure  1.  Example  of  Rademacher  functions 

These  functions  are  easy  to  generate  and  are  orthogonal, 
but  do  not  form  a  complete  basis.  They  therefore  cannot  be 
used  to  represent  a  general  y(t). 

The  orthogonality  can  be  shown  by  considering  any  two 
distinct  Rademacher  functions  R(n,t)  and  R(m,t),  where  n  and  m 
are  sequence  numbers  of  the  Rademacher  functions.  By  the 
construction  of  the  series,  the  number  of  intervals  for  each 
function  is  a  power  of  two.  Each  interval  of  constant  value 
is  of  equal  length  for  a  given  function.  Let  m  be  greater 
than  n.  Then  each  interval  of  fixed  value  for  R(n,t)  will 
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correspond  to  an  avan  number  of  equal  length  intervale  with 
alternating  eigne  for  R(m,t),  eo  the  integral  of  the  product 
of  the  two  functions  on  each  interval  will  be  eero. 

The  incompleteness  of  a  Rademacher  basis  can  be  seen  by 
noticing  that  there  is  always  an  odd  number  of  sign  changes. 
Using  Rademacher  functions  corresponds  to  restricting  Fourier 
analysis  solely  to  the  use  of  unshifted  sine  functions.  In 
fact,  one  way  of  defining  Rademacher  functions  is 

R(n,t)  ■=  sgn[  sin(  2nnot  )  ], 
where  n  is  the  sequence  number  of  the  Rademacher  function 

sgn[  *  ]  -  {  *_l  ;  ?  0. 

Any  cosine  function  of  the  form  cos(2not)  is  orthogonal  to  the 
basis.  Attempting  to  complete  the  basis  by  supplementing  it 
with  terms  sgn[  sin(  n2not  )  ],  where  n  is  not  a  power  of  2, 
and  sgn[  cos(  m2nut  )  ],  where  m  is  any  positive  integer, 
destroys  the  orthogonal  property. 

Since  Rademacher  functions  do  not  have  the  requisite 
properties  of  completeness  and  orthogonality,  they  cannot  be 
used  as  a  basis  for  spectral  analysis. 

2 . 3  Walsh  functions 

Walsh  functions  are  a  set  of  discrete-valued  functions 
which  assume  only  the  values  {-1,+1}.  They  are  orthogonal  and 
complete,  having  both  even  and  odd  symmetry.  They  can  be 
constructed  recursively,  as  products  of  Rademacher  functions, 
or  by  constructing  Hadamard  matrices  and  sorting.  Details  of 
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construction  will  b»  presented  later. 

Each  function  la  defined  by  convention  on  a  fixed 
interval  t  c  [0,T],  and  ia  written  WAL(n,t) .  As  with  the 
Rademacher  functions  T  is  usually  assumed  to  be  1  but  can  be 
any  value  if  the  function  is  scaled  appropriately.  The  value 
n  is  an  index  which  corresponds  uniquely  to  the  average  number 
of  zero  crossings  on  the  interval,  which  i6  called  the 
seguency  of  the  Walsh  function  in  an  analogy  to  frequency  in 
trigonometric  functions.  Walsh  functions  are  paired  by  even 
and  odd  symmetry  and  referred  to  as  CAL  and  SAL  functions, 
respectively.  These  are  defined  as  follows: 

CAL(k,t)  -  WAL(2k,t)  1  k  -  _ 

SAL(k,t)  *  WAL(2k-l,t)  J  * 

where  k  is  the  6equency.  Figure  2  illustrates  four  Walsh 
functions.  The  fifth  and  sixth  Walsh  functions  are  shown  to 
illustrate  that  Walsh  functions  are  not  periodic  in  their 
variations. 
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Walsh  and  Fouriar  spectra  havs  a  one  to  one 
correspondence.  This  can  be  seen  from  the  following  sets  of 
relations,  which  are  used  to  represent  a  tine  series  y(t)  in 
terns  of  a  Fourier  and  a  Walsh  basis,  respectively: 

CD  QD 

(Fourier)  y(t)  *  aQ  4  ][  a^cos  2nmt  +  £  bBsin  2nmt 

m«l  ®*1 

CD  CD 

(Walsh)  y (t)  -  Afi  +  X  A  CAL ( n , t )  4  £  BnSAL(n,t) 

n«l  n«l 

where  a^  bi#  A^ ,  and  are  the  Fourier  and  Walsh 
coefficients,  respectively,  denoted  as  cn'6  in  section  2.1. 
From  this  set  of  equalities  it  is  shown  by  Maqusi^11^  that 


cuiu 

®  p  1 

A  *  I  a-  cos  2nmt  CAL(n,t)  dt 
m=l  J0 

®  p  1 

B„  ■=  V  t  sin  2nmt  SAL(n,t)  dt 

or 

®  -1 

affi  *=  £  An  j  CAL(n,t)  cos  2nnt  dt 

®  rl 

bm  -  I  Bn  J  SAL(n, t)  sin  2nnt  dt 


CAL(n,t)  cos  2nmt  dt 


n  *  1,® 


SAL(n,t)  sin  2nmt  dt 


m  *  1,®. 


Clearly  Walsh  and  Fourier  spectra  can  be  converted  back 
and  forth.  We  have  found  the  linear  change  of  basis 
transform,  regardless  of  the  fact  that  both  bases  are 
infinite.  From  linear  algebra  we  know  that  for  any  complete 
series  we  consider  we  will  be  able  to  find  a  transformation  to 
the  familiar  Fourier  spectrum.  By  considering  alternative 
bases  we  have  neither  gained  nor  lost  information  relative  to 
a  Fourier  based  spectral  analysis.  However,  there  are  gains 
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to  be  had  by  using  tha  Walsh  rsprsssntation. 

The  first  advantage  for  Walsh  analysis  lies  in  the  number 
of  terms  needed  to  approximate  a  discrete  or  discontinuous 
time  series  from  the  simulation  output.  Walsh  functions  are 
more  efficient  than  trigonometric  functions  for  this  purpose, 
since  it  generally  requires  a  large  number  of  trigonometric 
terms  to  adequately  approximate  a  discontinuity.  For  the 
purposes  of  altering  the  value  of  a  discrete  input  parameter 
of  the  model,  trigonometric  functions  are  inappropriate. 

The  second  advantage  is  computational.  According  to 
Beauchamp  it  takes  n  log2n  COttPlex  mult ini icat ions  and 
additions  to  evaluate  a  fast  Fourier  transform.  Since  the 
Walsh  function  can  only  have  the  values  +l  or  -1,  evaluation 
requires  only  n  log2n  additions  if  the  analogous  fast 
transformation  is  used.  Although  the  degree  of  improvement  in 
efficiency  i6  machine  dependent,  addition  is  much  faster  than 
multiplication  on  digital  computers  using  current  technology. 


2.4  Haar  functions 


Haar  functions  are  discrete-valued  functions  which  are 


defined  as 


HAR(0,t)  =  1 
HAR(l.t)  =  { 

HAR(2P+  n,t)  - 


V  2P 
-V  2P 


0  £  t  £  1 

0  £  t  <  1/2 
1/2  £  t  £  1 

n/2p  £  t  <  (n  +  1/2) /2P 

(n  +  l/2)/2p  £  t  <  (n  +  1)/2P 

elsewhere. 


n 


This  definition  Allows  for  s  sequential  unique  numbering 
system.  The  first  seven  Hear  functions  are  illustrated  in 
figure  three.  Notice  that  there  is  no  unique  sequency 
correspondence  such  as  exists  with  the  Walsh  functions. 


Figure  3.  Example  of  Haar  functions 

A  different  indexing  scheme  groups  Haar  functions  by 
degree  i,  the  number  of  zero  crossings  in  width  2-1.  The 
definition  then  becomes 


22 


HAR  ( 0 , 2  ,  t )  -  1  0  £  t  £  1 

V'  21  (j  -  l)/2i  S  t  S  (j  *  1/2J/21 

HAR(i,j,t)  -  ■  W  2P  (j  -  2/2)/21  <£  t  <  j/21 

0  elsewhere 

2  *  0,1,2,...  j  ■  1, . . .,2. 

Only  one  tern  of  each  degree  can  be  utilized  if  Haar 
functions  are  used  for  identification.  Otherwise  it  would  be 
inpossible  to  distinguish  whether  an  observed  outcone  was  due 
to  a  specific  paraneter  or  due  to  a  lag  effect  fron  a 
different  paraneter.  This  property  could  actually  be 
advantageous  for  the  purpose  of  identifying  lagged  nodels. 
Furthermore,  there  is  a  fast  Haar  transform  which  requires 
only  2  (N  -  1}  additions. 

These  gains  are  offset  by  two  disadvantages.  The  first 
i6  that  Haar  functions  consist  of  three  states  and  so  cannot 
be  used  as  input  for  binary  parameters.  The  second  and  major 
disadvantage  is  that  there  is  not  a  convenient  product 
relationship  for  Haar  functions  such  as  exists  for  Walsh  and 
trigonometric  functions.  The  product  of  two  Haar  functions  is 
either  zero  if  there  is  no  overlap  of  the  non-zero  intervals, 
or  +  A^  HAR(k,t),  where  A^  is  the  amplitude  of  the  Haar 
function  with  the  larger  non-zero  interval,  and  HAR(k, t)  is 
the  function  with  the  smaller  non-zero  interval.  This  makes 
the  use  of  Haar  spectra  to  identify  interaction  terms  or 
higher  order  polynomial  terms  infeasible. 

The  remainder  of  this  report  will  concentrate  on  Walsh 
functions  since  they  appear  to  have  the  most  desirable 
properties. 
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It  should  be  noted  that  both  Welsh  and  Hear  analyses  are 
computationally  stable  procedures,  since  they  involve  only 
addition  operations.  This  is  not  true  of  Fourier  analysis. 
Addition  is  order  preserving  for  finite  precision  numbers, 
while  multiplication  is  not. 


At  thi6  point  ve  must  make  the  transition  from  continuous 
Walsh  functions  into  the  discrete  domain.  This  is  done  by 
scaling  the  time  axis  relative  to  the  highest  order  Walsh 
function  of  interest,  and  sampling  the  value  of  the  continuous 
function  over  unit  intervals.  We  will  end  up  with  a  vector  of 
N  numbers  which  correspond  to  sampling  WAL(k,t)  at  N  equal 
intervals,  i.e.  at  spacings  of  T/N.  Thus  WAL(k,i)  =  WAL(k,t) 
where  i  i6  the  integer  portion  of  [(Nt/T)+1],  WAL(k, . )  will 
be  used  to  denote  the  vector  consisting  of  WAL(k,i)  for 
i  «  1, . . . ,N. 

The  value  of  N  must  be  chosen  so  that  each  Walsh  function 
has  a  unique  vector  associated  with  it.  If  k  is  the  largest 


sequence  number  we  wish  to  observe  then  we  set  N  *  2 


flog  kl 


where  the  log  is  base  2  and  fx]  is  the  smallest  integer 
greater  than  x.  As  an  example,  if  we  were  interested  in  Walsh 
functions  up  to  order  5,  we  could  represent  them  as  vectors  of 
length  8  whose  elements  assumed  the  values  ±1. 

WAL ( 0 , . )  -  (  1,  1,  1,  1,  1,  1,  1,  1) 

WAL ( 5 , . )  -  (  1,-1, -1,  1,-1,  1,  1,-1) 


a 

s 


Arithmetic  operation*  on  Waleh  function*  make  extensive 
u*e  of  the  •  operator,  which  i*  called  a  dyadic  sum.  The 
dyadic  sum  is  a  bitwiae  XOR  operation,  where  p  XOR  q  is 
defined  by  the  following  table. 


examples  of  9  operator: 

7  111  9  1001 

•  5  <*>  C  101  9  3  <«>  9  0011 

2  010  10  1010 

We  present  a  short  list  of  interesting  properties  of 

Walsh  functions.  The  list  is  not  comprehensive.  The 

following  definitions  will  be  used. 

N  is  the  number  of  observations 

is  the  ith  parameter  in  the  model 
Yj  is  the  jth  term  in  the  output  time  series 

Y(X)  is  the  Xth  term  in  the  output  series  after  transformation 
Y  (X)  and  Y  (X)  are  the  X^  CAL  and  SAL  terms,  respectively 

C  D 

P(X)  is  the  spectrum  estimate  of  power  at  6equency  X 

It  is  worthwhile  to  note  here  that  the  Y's  above  are  the 
estimates  for  the  A's  and  B's  used  in  the  series 
representation  in  section  2.3.  We  are  using  the  notation 
y  h  Y  to  emphasize  that  the  original  observations  and  the 
estimated  coefficients  are  a  transform  pair,  either  of  which 
could  be  used  to  fully  reconstruct  the  other. 


■) 
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PROPERTIES 

symmetry 

multiplication 

properties 


Walsh 

transformation 


spectrum 

estimator 


WAL(k,i)  -  WAL(i,k) 

WAL ( k ,  i ) WAL ( j , i )  •  WAL ( j  ,  i ) WAL ( k ,  i ) 

WAL ( k ,  i ) WAL ( j  ,  i )  »  WAL(k  ©  j,i) 
CAL(k,i)CAL(j,i)  -  CAL(k  ©  j,i) 
SAL(k,i)CAL(j,i)  -  SAL( j  ©  [k-l](i) 

SAL ( k , i ) SAL ( j , i )  «  CAL( [k-1  ]  ©  [ j-1 ] , i) 
.  N-l 

Y(k>  "  s  iIoyiWAL(k»1) 


P(0)  *  Y*(0> 


P(k)  -  Y^(k)  +  Yg(k) 


k= 1  #  2  ,  .  .  .  ,  ( N/ 2 ) — 1 


y;(n/2> 


to  ■ 

h” 


P(H/2) 

There  is  a  shift  theorem  for  trigonometric  functions 
which  states  that  shifting  a  trig  function  does  not  change  the 
observed  frequency.  There  is  no  comparable  theorem  for  Walsh 
functions.  Their  shift  behavior  is  discussed  later. 

Two  immediate  results  for  the  ®  operator  are  that 
k  ©  k  =  0  for  all  k;  and  k  ©  0  ■  k  for  all  k.  Applying  these 
to  the  Walsh  multiplication  properties  shows  that  if 
x  =  WAL(k,i) ,  then 


w 

V,*  *\'  •*.*  *’  *  •{J  /•*./ 


O' 

W_*-  - 


- 


X  •  X  *  X  •  WAL(k  t  k,i)  •  WAL(0, 1) 

X3  -  X  *  X2  -  WAL(k  *0,1)  -  WAL(k,l) 

x4  -  X  *  x3  -  WAL(k  •  k,l)  -  WAL (0,1) 

• 

By  induction  all  even  powers  of  x  will  have  an  indicator 
sequency  of  0  and  all  odd  powers  will  have  sequency  k/2.  This 
Deans  that  we  will  be  unable  to  estimate  any  but  first  order 
and  interaction  terms  using  a  single  Walsh  function.  If  the 
problem  is  viewed  geometrically  this  is  reasonable,  since 
using  a  Walsh  function  corresponds  to  sampling  the  data  at 
only  two  points.  With  two  points  it  is  not  possible  to 
determine  more  than  a  first  order  model. 

If  we  are  only  interested  in  significant  parameter 
detection,  this  presents  no  difficulty.  However,  if  we  hope 
to  use  the  data  later  to  estimate  a  model  there  may  be 
trouble.  The  situation  is  acceptable  in  the  case  of  binary 
parameters,  but  there  is  a  prc^em  for  continuous  or  p-state 
parameters  where  p  >  2.  We  would  prefer  to  have  a  single 
procedure  which  works  for  both  continuous  and  discrete 
parameters.  One  way  to  overcome  the  problem  is  to  use  sums  of 
Walsh  functions  as  the  input.  A  sum  of  two  Walsh  functions 
with  amplitudes  A^  and  A2  will  yield  {+A^  +A2),  or  up  to  four 
distinct  sampling  points.  This  should  be  sufficient  to 
construct  most  non-linear  models  in  practice. 
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X ( i )  -  WAL(j,i)  4  WAL(k, 1) 
y(l)  -  h1  4  a2x2(i) 

-  4  a2[WAL(j,l)  4  WAL(k,i)]2 

-  a1  4  a2[WAL(j  9  j , i)  4  2WAL( j  9  k,i) 

4  WAL(k  9  k,  i)  ] 

-  (ax  4  2a2)WAL(0,i)  4  2a2WAL(j  9  k,i) 

We  would  observe  an  Increase  in  the  steady  state  term, 
sequency  0,  and  an  interaction  term  at  6equency  (j  9  k)/2  for 
the  square  term  in  a  model  of  this  form.  A  linear  term  in  x 
would  just  have  sequencies  j/2  and  k/2  appear  in  the  spectrum. 

This  form  of  input  may  offer  a  viable  approach  to 
detecting  system  gain.  Gain  is  the  system's  tendency  to 
amplify  or  attenuate  the  response  at  different  sequencies.  If 
a  parameter  has  a  linear  term,  in  the  absence  of  system  gain 
the  amplitudes  of  WAL(j,i)  and  WAL(k,i)  should  be  the  same.  By 
using  a  combination  of  high  sequencies  and  low  sequencies,  and 
comparing  the  spectra  of  the  two  sequencies  for  a  given  run, 
any  differences  in  power  should  be  due  to  system  gain  or 
stochastic  error.  We  can  statistically  test  the  spectrum  for 
gain  by  performing  an  F-test  as  described  later. 

According  to  the  literature  Walsh  functions  are 
theoretically  better  than  trigonometric  functions  for  parameter 
detection  in  non-linear  models.  This  can  be  explained  by 
figure  4.  Since  input  is  held  at  a  single  value  for  each  time 
interval,  the  system  output  has  the  characteristic  of 
maintaining  a  single  value  over  the  same  interval.  The 


amplitude  of  tha  output  la  changed,  but  the  number  of  ateps 
remalne  the  eame  and  will  be  detected  by  Walsh  analysis. 


Otit  put 


Figure  4.  A  nonlinear  response  surface 

A  problem  in  the  use  of  Walsh  analysis  is  that  it  is  not 
invariant  under  phase  shifts.  According  to  Beauchamp  the  CAL 
and  SAL  functions  of  a  given  seguency  vary  inversely  to  each 
other  in  such  a  manner  that  the  spectrum  estimator  is 
"relatively  insensitive"  to  changes  in  phase.  However,  the 
practitioner  should  be  cautious  if  the  simulation  model 
involves  time  lags.  Our  experiments  indicate  that  the  power 
spectrum  is  phase  shift  invariant  if  the  seguency  is  N/2  or 
N/4,  where  N  is  the  number  of  terms  in  the  series.  For  other 
seguencies  the  spectrum  will  have  a  spike  at  the  driving 
seguency,  but  there  will  be  spikes  of  varying  heights  at  other 
seguencies  as  well.  The  total  power  in  the  spectrum  remains 
constant,  and  empirical  evidence  indicates  that  the  bulk  of 
the  power  is  usually  displayed  at  the  original  seguency. 
Spectra  are  usually  plotted  on  a  log  scale,  and  it  is 


difficult  to  discsrn  the  primary  from  ••condary  spikes. 
Furthermore,  thara  ara  caaes  whara  tha  power  ia  uniformly 
distributed  across  all  spikes,  making  identification 
impossible.  We  have  been  unable  to  find  an  analytical  method 
of  predicting  the  location  of  the  extra  spikes,  other  than  by 
actually  producing  spectra  of  the  aeguency  at  a  specified 
lag.  See  appendix  2  for  details. 

2 . 6  Generating  Walsh  functions 

The  vector  notation  introduced  in  the  previous  section 
allows  us  to  use  matrices  to  represent  Walsh  functions,  which 
simplifies  notation  substantially .  Throughout  the  remainder 
of  this  paper  we  will  use  WN  to  denote  an  NxN  matrix  comprised 
of  the  first  N  Walsh  functions. 
example 


1  1  1  1  1 


L  1  -1  1  -1  J 


Walsh  functions  can  be  generated  recursively  by  the  use 
of  Hadamard  matrices.  One  definition  of  a  Hadamard  matrix  is 
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H2(n4l) 


H 


H 

2 


n 


J 


(n  £  0)  . 


example 


1  1  1  1  1 


Notice  that  Hv,  contains  the  first  N  Walsh  functions  as  its 
rows  and  columns,  but  not  in  seguency  order.  After  sorting  we 
have  W,  which  retains  the  symmetry  property  of  H.  However, 
while  thi6  i6  an  interesting  mathematical  viewpoint  it  does 
not  provide  a  convenient  and  efficient  method  for  obtaining 
Walsh  functions. 

The  discrete  Walsh  transformation  can  be  viewed  as  a 
matrix  transformation  using  W  which  performs  a  change  of  basis 
on  a  vector  y  of  length  N  =  2n. 

Y  .  i  H  y 

By  the  orthonormal  property 

W  WT  -  N  I 

where 

T  is  the  transpose  operator 
N  is  as  defined  above 
I  is  the  NxN  identity  matrix 

T 

W  =  W  since  W  is  symmetric.  W  is  invertible  because  it  is 
composed  of  orthogonal  vectors,  and  so  must  be  of  full  rank. 
We  can  compute  W  1  as  follows: 


21 


-1  T  -1 

W  AW  W  »  N  W  J 

T  -1 

W  «  N  W  a 

-1  IT 

«  -  8  M 

K*1  -  i  W. 

This  means  that  the  Walsh  transformation  is  its  own 
inverse  after  scaling  by  N.  This  provides  us  with  perhaps  the 
simplest  method  of  generating  Walsh  functions.  Just  as  with 
the  Fourier  transformation  matrix,  W  has  redundancies  and  can 
be  reduced  to  a  sparser  matrix.  Use  of  this  sparse  matrix 
results  in  a  considerable  savings  in  computational  time  and 
complexity.  This  is  called  the  Fast  Walsh  Transform  (FWT) . 

By  performing  a  FWT  on  a  vector  which  has  the  value  one  in  the 
location  corresponding  to  the  desired  sequence  number  and 
which  is  zero  elsewhere,  we  obtain  the  indicated  Walsh 
function  as  output, 
example 


WAL ( 0 , .) 


=  WAL ( 2 , . ) 


A  computer  algorithm  for  the  FWT  is  included  in  appendix  1, 


Section  3 


Dtelgn  of  a  Spectral  Experiment 


3.1  DlgtrlfrvtUpn  -pi  .the,  fipsctrum  Jptlflftfcgr 

We  assume  that  the  output  from  a  simulation  of  a 
stochastic  system  is  composed  of  a  deterministic  signal  term  s 
and  an  error  term  c . 

*i  -  »i  +  ei 

We  will  show  that  the  spectrum  estimator  has  a  x  distribution 

under  the  assumption  that  the  additive  error  term  has  the 

distribution  of  discrete  white  noise  as  defined  in  Jenkins  and 

Watts.  This  means  that  it  i6  normally  distributed  with  mean 

2 

zero,  variance  o£  ,  and  cov(c^,Cj)  «  0  for  i  f  j .  Using  the 
terminology  defined  in  section  2.5  we  obtain 

Y«|wy«|w  (s  +  c  )  «£ws  +  £we. 

The  random  component  associated  with  this  is 


Since  the  only  variability  in  Y  (i)  is  due  to  c.,  clearly  the 

variance  of  ¥, (i)  is  of  and  Cov(Y„ (i) , ¥_  ( j ) )  *  0  for  i  *  j . 
c  c  c  c 

Since  Y  (i)  is  a  linear  function  of  a  normal  random 
variable,  it  is  distributed  normally.  Thus 

Yf(i)  Yf(k) 

Var(Yc(i))  of 

2 

has  a  X  distribution  with  v,  the  degrees  of  freedom,  equal  to 

one.  Since  the  covariance  of  the  Y  's  is  zero  the  Y  (i)'s  are 

c  c 

independent . 


This  enables  ua  to  parfom  a  statistical  taat  of  whether 
two  independent  spectrum  estimates  ere  equal.  Hence  we  can 
test  whether  an  observed  spike  in  the  spectrum  is  due  to  a 
factor  effect  or  is  due  to  error. 

The  next  issue  is  how  to  increase  the  degrees  of  freedom 
for  the  estimator.  Consider  a  run  in  which  all  of  the 
sequencies  assigned  are  less  than  or  equal  to  k.  Let 

K  «  2  f  where  the  log  is  base  2.  Then  all  crossterms 

will  also  be  of  order  less  than  K,  since  the  •  operator  works 
in  a  bitwise  manner.  We  therefore  can  restrict  the  analysis  to 
batches  of  K  points.  If  we  now  let  the  simulation  run  for  m 
batches  of  length  K  we  will  have  m  independent  estimates  of 
each  of  the  K  points  under  the  assumption  that  the  error  is  an 
additive  discrete  white  noise  term,  since  we  are  performing 
linear  transformations  on  independent  observations.  We  can 
now  construct  K  new  estimators  from  the  sum  of  these  m 

2 

estimates.  The  new  estimates  will  be  distributed  as  x  random 

variables  with  v  ■  2m  except  at  the  endpoint6,  since  they  are 

2 

constructed  as  a  sum  of  m  x  variables  with  v  *  2.  The 

2 

endpoint  estimators  are  a  sum  of  m  x  distributions  with 
v  *  1,  and  hence  will  have  v  «  m. 

3 . 2  Seauencv  Selection 

The  basic  problem  is  to  find  a  set  of  sequencies  for  n 
factors  such  that  all  of  the  original  terms  and  the  crossterms 
have  a  unique  sequency.  A  straightforward  solution  is  to 
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assign  the  factor#  eequenciee  which  are  power#  of  two,  which 
actually  assigne  one  bit  to  aach  factor.  Thie  will  always 
yield  a  unique  eet  of  croee  terms,  and  would  be  ideal  if  we 
were  interested  in  all  possible  interactions.  For  example, 
for  a  three  factor  model  with  factor  i  assigned  sequency  2i-1 
we  might  observe 

sequency  -»  term 


1 

2 

4 
3 

5 

6 
7 


If  we  are  interested  only  in  two-term  interactions  this 

is  inefficient.  For  a  first  order  model  there  are  n  original 

2 

terms  and  n  choose  2  crossterms,  or  a  total  of  (n  +  n)/2 

sequencies.  We  want  to  minimize  the  required  number  of 

v+i 

observations,  and  the  FWT  needs  2  observations  to  detect  a 

•  V — 1  V  V+ 1 

sequency  p  in  the  range  2  s  p  <  2  .  The  number  is  2 
rather  than  2  because  two  Walsh  functions  are  used  for  each 


sequency  estimate.  We  can  obtain  substantial  savings  even  for 

2 

relatively  small  n  if  we  can  find  (n  +  n)/2  terms  rather  than 

2n  terms.  For  example,  if  we  have  10  factors  we  will  generate 

55  sequencies.  If  we  select  the  original  10  in  such  a  way 

that  all  55  sequencies  are  less  than  64,  then  we  need  only  128 

observations.  This  is  much  more  efficient  than  taking 
11 

2  (*2048)  observations  to  identify  the  same  model. 


Recall  from  section  2.6  that  in  order  to  detect  a  square 
tern  ve  need  two  Walsh  seguenciss.  If  the  square  tern  is 
present  we  will  observe  the  crosstem  of  the  two  seguencies  in 
the  output.  We  can  solve  the  selection  problen  as  if  we  had  a 
nodel  with  twice  as  nany  factors,  and  assign  two  unique 
sequencies  per  factor.  Alternatively  we  can  choose  the 
sequencies  as  if  there  were  one  extra  factor  and  nake  one 
sequency  connon  to  all  factors  being  tested  for  a  square 
tern.  Each  of  the  two  proposals  has  nerit.  The  first  one 
allows  a  pairing  of  high  and  low  sequencies  for  every  factor 
so  that  the  nodel  can  be  tested  for  systen  gain.  The  second 
has  the  desirable  property  of  reducing  the  size  of  the  problem 
ve  must  solve  to  assign  sequencies.  The  user  must  decide 
which  of  these  considerations  is  more  important  for  his 
specific  application. 

In  order  to  solve  the  proposed  problem  we  have  written  a 
computer  program  called  SEQ  which  enumerates  sets  of  n 
sequencies,  generates  the  n  choose  2  crossterms,  sorts,  and 
checks  for  uniqueness  of  all  elements.  It  prints  out  the 
original  n  sequencies  if  they  constitute  a  solution.  The  code 
can  be  found  in  appendix  1. 

If  we  did  this  by  explicit  enumeration  of  all 
configurations  the  problem  would  rapidly  become  too  large  to 
solve.  To  find  3  input  sequencies,  for  example,  we  have  6 
terms.  We  want  these  to  fall  in  the  range  1  to  7.  We  can 
without  loss  of  generality  order  them  so  that  the  first  is 
smallest  and  the  third  is  largest.  Then  to  have  uniqueness  the 
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first  can  only  rang*  from  1  to  5,  tha  aacond  from  2  to  6,  and 
the  third  from  3  to  7 .  If  we  call  tha  number  of  parameters  n 
and  the  highest  sequency  of  interest  is  m,  in  our  example 
n  ■  3  and  m  ■  7.  Then  the  number  of  configurations  possible 
is  m  choose  n.  For  our  example  this  is  35.  For  n  ■  10  we 
need  at  least  (n  +  n)/2  ■  55  seguencies.  We  might  as  well 
consider  all  sequencies  up  to  the  next  power  of  two,  i.e. 
seguencies  0  through  63,  since  they  will  all  yield  the  same 
batch  size.  Then  there  are  63  choose  10  configurations,  or 
more  than  1011  possibilities  to  consider.  While  there  may  be 
several  solutions,  it  could  take  a  long  time  to  find  any  of 
them. 

We  can  drastically  reduce  the  amount  of  work  required  by 
noting  that  for  the  set  {x1>*2 , . . .  ,xn)  to  be  a  solution  the 
subset  {*ltx2 , . . .  »xn_1)  must  also  be  a  solution.  Hence  if  we 
can  find  a  solution  to  the  problem  at  level  k,  we  can  add  one 
additional  element  and  range  through  the  possible  values  until 
we  have  a  solution  at  level  k+1.  If  we  fail  to  find  any 
solutions  we  return  to  level  k  and  search  for  a  new  solution 
there.  By  proceeding  iteratively  in  this  fashion  we  can  build 
the  solutions  to  fairly  large  problems  quite  easily.  The 
program  presented  in  appendix  1  can  handle  the  case  n  ■  16  in 
under  three  minutes  on  an  IBM  personal  computer.  Some 
solutions  for  various  values  of  n  are  presented  here. 
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n  • 

4 

{ 1 , 2 , 4 , 8  ) 

n  • 

5 

{1,2,4,8,15) 

n  - 

6 

{1,2,4,8,15,16) 

n  ■ 

16 

{1,2,4,8,15,16,32, 
51,64,85, 106,128, 
150,171,219,237) 

A  separate  program  called  XTERM  is  included  for 
predicting  all  the  crossterms  from  a  given  set  of  input 
6eguencies .  Sample  output  from  that  program  is  shown  here  for 
the  case  n  ■  4. 


enter  the  number  of  seguencies: 
enter  the  sequencies: 


1  xor  2:  3 

1  xor  4:  5 

1  xor  8:  9 

2  xor  4:  6 

2  xor  8:  10 

4  xor  8:  12 


8 


4 


3 . 3  Design  of  an  Experiment 

We  will  now  describe  a  general  methodology  for  designing 
a  Walsh  spectral  experiment  for  identification  purposes. 

The  experimenter  must  first  establish  the  number  of 
factors  to  be  identified,  and  the  order  of  the  model  to  be 
fit.  If  a  model  of  greater  than  order  2  is  desired  we 
recommend  using  sequencies  which  are  powers  of  two,  as 
described  in  section  3.2.  Otherwise  the  experimenter  can  use 
the  programs  included  in  appendix  1. 
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Next  the  experimenter  should  decide  whether  he  needs  two 
■equencies  per  factor  or  can  utilise  a  common  term  approach. 
Two  sequencies  per  factor  rapidly  increases  the  required 
number  of  sequencies,  and  hence  increases  the  number  of 
observations.  For  a  minimum  number  of  observations  the  common 
term  approach  is  preferable,  but  it  has  two  disadvantages. 

The  first  is  that  it  may  be  more  susceptible  to  system  gain 
behavior,  since  not  every  factor  can  be  assigned  both  a  high 
and  a  low  sequency.  The  second  disadvantage  is  that  it  is  not 
easy  to  discern  the  difference  between  crossterms  and  square 
terms  in  the  model.  For  example,  consider  the  following  two 
factors 

xa  -  a2  [WAL(p,t)  +  WAL ( q , t ) ] 
x2  -  a2  [WAL(p, t)  +  WAL(r,t)] 

where  aj^  and  a2  are  the  amplitudes  of  x^^  and  x2  respectively. 
Then 

X1  "  al  [WAL(P«P/t)  +  WAL(q©q,t)  +  2  WAL(p©q,t)] 

«  2  a*  [WAL(0,t)  +  WAL(p©q,t)] 

X1X2  “  a^lWAL^t)  +  WAL(q,t)]  [WAL(p,t)  +  WAL(r,t)] 

*  a^JWALfpep,^  +  WAL(p©r,t)  +  WAL(p©q,t)  +  WAL(q©r,t)] 

«  a1a2[WAL(0,t)  +  WAL(p©r,t)  +  WAL(p©q,t)  +  WAL(q©r,t)] 

2 

Both  the  x1  and  x^2  terms  share  WAL(p©q,t)  as  a  common 
identifying  sequency.  In  the  absence  of  terms  WAL(pQr)  and 
WAL(q©r)  we  would  conclude  that  the  spike  in  the  spectrum  was 
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due  to  an  affact,  but  if  the  other  terse  are  present  we 

cannot  exclude  the  poseibility  that  there  la  an  x^  affect  as 

veil  as  an  x^x2  effect.  Note,  however,  that  all  three  spikes 

identifying  an  x,^  tern  should  have  the  sane  height.  If  the 

height  at  p©q  ie  greater  than  the  other  two  heights,  it  is 

2 

probably  due  to  the  addition  of  an  x^^  tern,  although  it  night 
be  due  to  gain. 

After  deciding  on  the  number  of  aeguencies  which  will  be 
needed,  the  experimenter  should  select  a  set  such  that 
crossterms  are  uniquely  identifiable.  The  program  SEQ  in 
appendix  1  can  be  used  for  this  purpose. 

The  user  should  then  run  the  experiment  by  varying  each 
factor  according  to  the  even  numbered  Walsh  function  which 
corresponds  to  its  assigned  sequency.  The  odd  sequence  number 
cannot  be  used  because  it  does  not  preserve  the  ©  operator 
under  which  the  sequencies  were  selected,  while  the  even 
sequence  numbers  do.  For  example,  if  parameters  x1  and  x2  are 
assigned  sequencies  4  and  8,  respectively,  then  we  have  the 
following. 

Seq(4)  <=>  WAL(7,t)  ,  WAL(8,t) 

Seq(8)  <«>  WAL(15, t)  ,  WAL(16,t) 

If  the  even  terms  are  used 

WAL(8,t)  WAL ( 1 6 , t )  *  WAL(8©16,t) 

-  WAL ( 2  4 , t ) 

*>  seq(12) 

*  seq(8©4) . 
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If  the  odd  term*  are  used 

WAL(7,t)  WAL( 15 , t )  -  WAL(7§15,t) 

-  WAL(8 , t ) 

->  aeq(4) 
x  seq(8©4) . 

Thus  ve  can  only  use  the  sequency  based  experimental  design  if 
we  convert  to  even  numbered  Walsh  functions. 

Finally,  the  output  from  the  experiment  should  be 
evaluated  using  a  program  such  as  SPECTRUM  in  appendix  1.  The 
spectrum  obtained  from  that  program  can  be  analyzed  for 
sequency  components  corresponding  to  the  factor  inputs. 
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Section  4  Analysis  and  Concluaiona 

4.1  Veriflcatlon_gf  Walsh  Analysis 

Before  proceeding  further  we  wanted  to  verify  the  ability 
of  Walsh  analysis  to  extract  the  correct  signal  in  the 
presence  of  noi6e.  To  this  end  we  constructed  a  file  of 
random  numbers  to  be  used  as  a  signal  with  additive  noise.  A 
second  file  consisted  of  the  fractional  parts  of  the  numbers 
from  the  first  file,  and  was  considered  as  the  noise.  A  third 
file  was  constructed  from  the  integer  portion  of  the  first 
file,  and  considered  to  be  a  noiseless  signal.  The  three 
files  are  presented  here. 

Signal 

4.1  3.9  4.1  3.9  3.8  4.7  5.2  6.4 

4.1  3.9  4.1  3.9  3.8  4.7  5.2  6.4 


Noise 

0.1  0.9  0.1  0.9  0.8  0.7  0.2  0.4 
0.1  0.9  0.1  0.9  0.8  0.7  0.2  0.4 


Noiseless 

4.0  3.0  4.0  3.0  3.0  4.0  5.0  6.0 
4.0  3.0  4.0  3.0  3.0  4.0  5.0  6.0 

All  three  files  were  subjected  to  spectral  analysis  using 
the  program  SPECTRUM  in  appendix  1.  Each  file  was  plotted 
individually,  and  the  signal  and  noise  files  were  used  to  form 
the  signal  to  noise  ratio.  The  spectral  plots  are  presented 
as  figures  5  through  8.  The  signal  plot  in  figure  5  has  a 
high  seguency  component  which  is  not  present  in  the  noiseless 
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plot  of  figur#  7.  This  noiss  component  le  successfully 
eliminated  by  taking  the  signal  to  noise  ratio  in  figure  6. 
That  plot  produces  spikes  in  the  same  location  aB  those  in  the 
noiseless  spectrum,  indicating  that  we  have  identified  the 
correct  seguency  components  of  the  true  signal  in  the  presence 
of  noise. 

The  noise  component  in  this  case  does  not  meet  the 
assumption  of  being  normally  distributed  with  mean  zero.  This 
indicates  that  the  discriminating  ability  of  Walsh  spectra  is 
not  necessarily  dependent  on  the  distribution  of  the  error 
terms.  That  distribution  is  important  only  in  being  able  to 
construct  a  statistical  test  for  identifying  signficant 
terms.  If  the  assumptions  of  section  3.1  hold,  the  values  to 
the  left  of  the  signal  to  noise  ratio  plot  are  F  values,  and 
can  be  used  directly  after  computing  the  appropriate  degrees 


of  freedom. 


enter  the  name  of  the  eignal  file: 
signal 

enter  the  name  of  the  noise  file  or  the  word  "nul": 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  16 


2 . 0362650E+01  0+** *************************************  * 

0. OOOOOOOE+OO  + 

4 .1281260E-01  +****************** 

0. OOOOOOOE+OO  -j- 

1.5156250E-01  +************ 

0. OOOOOOOE+OO  5+ 

9 . 9062  500E-02  +********** 

0. OOOOOOOE+OO  + 

4 . 5156280E-02  +***** 

Figure  5.  Spectrum  of  signal  run  with  noise  component 


enter  the  name  of  the  signal  file: 
noise 

enter  the  name  of  the  noise  file  or  the  word  "nul": 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  16 


2.6265620E-01  Of**************************************** 

0. OOOOOOOE+OO  f 

1.2812500E-02  -}-********** 

0. OOOOOOOE+OO  f 

1 . 4062500E-02  -j-*********** 

0. OOOOOOOE+OO  5+ 

3 . 6562490E-02  f ******************** 

0. OOOOOOOE+OO  + 

4 . 5156240E-02  f ********************** 


Figure  6.  Spectrum  of  noise  component 
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enter  the  name  of  the  signal  fila: 

noiseles 

enter  the  name  of  the  noise  file  or  the  word  *'nul": 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  16 


1 . 6000000E+01  o+* ***********************  **************** 

O.OOOOOOOE+OO  + 

5. 0000000E-01  -j-************* 

O.OOOOOOOE+OO  + 

2.5000000E-01  +******* 

O.OOOOOOOE+OO  5f 

2 . 5000000E-01  +******* 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 


Figure  7.  Spectrum  of  signal  without  noise  component 


enter  the  name  of  the  signal  file: 
signal 

enter  the  name  of  the  noise  file  or  the  word  Mnul” : 
noise 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  16 


7. 752587 0E+01  Of**************************************** 

1 .  OOOOOOOE+OO  4-******* 

3 . 22 19520E+01  -^*******  ******  ******************** 

1. OOOOOOOE+OO  4******* 

1 . 0777780E+01  +************* ************ 

1. OOOOOOOE+OO  5-j-******* 

2 . 7094020E+00  +*************** 

1. OOOOOOOE+OO  -}-******* 

1 . 0000010E+00  +******* 

Figure  8.  Spectrum  of  signal/noise  ratio 
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4.2  A  Polynomial  Model 

We  next  tested  a  model  with  4  factor*.  The  model  16  as 
follows: 


*t  “  xlft  +  xl,t  4 


+ 


X2,tX3,t 

3 


+  .001  x.  4.  +  error 
4 ,  t  t 


where  x^^  through  x4  are  factors 

errort  -  .8  errort_1  +  .6  z 
z  -  N (0 , 1)  . 

The  spectral  plot  included  as  figure  9  is  based  on  3  batches 
of  32  observations.  Sequencies  were  assigned  as  follows, 
using  the  common  term  scheme. 


Ll* 

V 


seq(15)  +  seq(8) 
seq(15)  +  seq(4) 


x3 :  seq(15)  +  seq(2) 

x4 :  seq(15)  +  seq(l) 

We  would  predict  the  following  sequencies  to  be  observed. 
xx  :  seq(8) ,seq(15) 
x,2  :  seq(0) ,seq(7) 
x1  :  seq(4) ,seq(15) 
x2x3:  6eq(0) ,seq(6) ,seq(ll) ,seq(13) 

x4  :  seq(l) , seq(15) 

The  derivations  of  the  linear  terms  should  be  evident.  We 
will  derive  the  x2x3  crossterm  as  an  example. 


x2x3  -  (••q(15)  4  ••q(4) j  [ssg(15)  4  saq(2)) 

-  (••q(15)2]  4  [ seq (15)  ssq(2)]  4  [itq(15)  »eq(4)] 

4  [ssq(4)  seq (2 ) ] 

■  seq(15©15)  4  seq(15«2)  4  seq(15©4)  4  seq(4©2) 

*  seq(O)  4  seq(13)  4  seq(ll)  4  seq(6) 

Spikes  were  observed  at  all  predicted  locations,  and  as 
expected  the  spike  for  x^  was  negligable.  The  spike  for  the 
Xj/3  term  was  small  but  noticable,  and  would  probably  be 
accepted  as  significant  if  more  observations  were  taken 


enter  the  nan#  of  tha  signal  fila: 

Big. run 

enter  tha  name  of  tha  noiaa  fila  or  tha  word  Mnul": 
naa . run 

antar  the  number  of  batches  and  batch  size  of  the  files: 

3  32 


4 .4337260E+01 
9 . 98  3  5490E-01 
9 . 9999980E-01 
9 . 9999940E-01 
1 . 502 1450E+00 
1 . OOOOOOOE+OO 
1.1210890E+01 
2.1937100E+02 
4 . 2658360E+01 
1. OOOOOOOE+OO 
1 . 000001 0E+00 
4.9413060E+01 
9 . 9999920E-01 
1 . 0608800E+01 
9 . 9999960E-01 
1 . 4 179650E+02 
9 . 9999980E-01 


Of****************************** 
f ****** 
f ****** 
f ****** 
f ******** 

5f ****** 

+****n********************************** 

f ****** 

10f ****** 

f ****** 

-j-*  ******************** 

-j-  ****** 

f****** 


Figure  9.  Spectrum  of  signal/noise  ratio  for  polynomial  model 
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Section  5 


We  present  here  several  topics  which  are  interesting 


extensions  of  the  current  work. 


IrSiTTiC 


LtJBnitMiIwMil: 


While  Walsh  and  traditional  spectral  analysis  can  be 
applied  to  any  time  series,  it  is  preferable  to  tailor  the 
analysis  to  the  type  of  system  being  analyzed.  One  area  for 
future  research  is  to  investigate  the  use  of  generalized  Walsh 
functions  for  analyzing  p-valued  parameters.  Elementary  Walsh 
functions  are  binary,  and  can  be  used  to  evaluate  binary 
parameters.  Generalized  Walsh  functions  are  series  which  have 
the  orthogonality  and  completeness  properties  while  assuming  p 
discrete  values,  where  p  £  3.  Using  these  functions  as  the 
basis  for  spectral  analysis,  it  would  be  theoretically 
possible  to  evaluate  systems  with  multi-valued  discrete 
parameters  by  spectral  methods. 


5.2  Parameter  Estimation 


Once  the  significant  factors  have  been  identified  a 
logical  next  step  is  to  attempt  parameter  estimation  for  our 
model  of  the  system.  This  problem  will  potentially  be 
complicated  by  non-linearities  in  the  system.  After  some 
preliminary  readings  we  would  recommend  an  iterative  least 
squares  approach. 
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On*  *xt*n*ion  of  parameter  estimation  would  b*  to  attempt 


to  reduce  th*  matrix  of  coefficienta  to  find  a  minimal 
representation  for  th#  I/O  behavior  of  the  aystem.  A 
representation  of  this  sort  would  have  many  possible 
applications.  It  could  be  used  for  atability/sensitivity 
analysis  of  the  system  via  analytic  methods.  It  could  also  be 
used  as  an  external  control  aystem  for  variance  reduction  of 
the  computer  simulation.  It  might  be  fruitful  to  investigate 
whether  there  i6  a  relationship  between  minimal  I/O  systems 
and  minimal  representations  of  the  simulation  using  event 
graphs . 

5.4  Parameter  Optimization 

Optimization  using  computer  simulation  i6  generally 
regarded  as  a  poor  idea.  It  quickly  becomes  an  expensive 
proposition  if  one  employs  the  traditional  approach  of  using 
each  run  a6  a  single  observation.  However ,  using  spectral 
methods  it  should  be  possible  to  use  simulation  for 
optimization  at  a  substantially  reduced  cost.  One  possible 
approach  would  be  to  try  to  converge  iteratively  within  a 
single  run.  Another  approach  would  be  to  try  to  solve  the 
problem  analytically  after  estimating  the  system. 
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of  the  system.  However,  It  1b  still  possible  to  improve 
variance  reduction  without  estimating  a  model  by  using  term 
sensitivity  to  indicate  which  variance  reduction  technique  is 
appropriate.  For  example,  antithetic  variates  should  only  be 
used  if  the  response  surface  is  monotonic,  so  the  presence  of 
non-linear  polynomial  terms  would  contraindicate  the  use  of 
this  method. 

5.6  Multiple-Input/Multlple-Qutput -Models 

Thus  far  all  work  done  by  Schruben,  Cogliano,  and  Sanchez 
has  focused  on  single  output  models.  Work  needs  to  be  done  so 
that  the  procedure  can  be  applied  to  more  general  simulation 
models. 

5.7  Statistical  Design  of  Experiments 

It  is  worth  noting  that  the  charts  found  in  many 

y 

experimental  design  texts  for  setting  factor  levels  in  a  2 
factorial  experimental  design  are  matrices  of  Walsh 
functions.  One  can  conjecture  that  a  spectral  approach  to 
experimental  design  has  been  overlooked  prior  to  now  because 
most  statisticians  are  familiar  only  with  the  traditional 
Fourier  spectrum.  It  may  be  possible  to  gain  additional 
insights  into  design  problems  by  using  discrete  spectra.  For 
instance,  generalized  Walsh  functions  may  be  useful  for 

y 

designing  and  analyzing  p  factorial  experiments. 
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Appendix  1 


All  program*  vara  implemented  In  tha  Paacal  language. 
Paacal  vaa  choaan  for  lta  readability  and  tranaportability . 
Standard  faaturaa  a r#  uaad  throughout  with  the  exception  of 
the  program  for  enumerating  aequencies.  That  program  makes 
use  of  the  IBM  and  Turbo  pascal  compilers  ability  to  do  XOR 
operations  on  integers.  Several  programs  were  developed  on 
the  Purdue  UNIX  system  and  then  implemented  on  an  IBM 
personal  computer. 

The  first  program  is  a  procedure  which  does  a  Fast  Walsh 
Transform  using  the  signal  flow  diagram  in  figure  3.3  of 
Beauchamp. 

The  second  program  generates  permutations  of  seguencies 
and  cross-terms  for  a  given  number  of  parameters.  After 
generation  the  sequencies  are  sorted  using  a  quicksort 
algorithm.  They  are  checked  for  redundant  values.  If  there 
are  no  redundancies  the  input  sequencies  are  printed. 

The  third  program  will  take  a  list  of  sequencies  as 
input,  and  generate  all  the  crossterms  as  output.  This 
procedure  is  described  in  section  3.2. 

The  fourth  program  is  a  general  procedure  for  doing 
Walsh  analysis  on  an  input  data  set.  It  prompts  the  user  for 
the  location  of  the  data,  number  of  batches,  and  size  of  each 
batch.  It  then  does  the  Walsh  analysis  and  prints  the  output 
as  a  barplot.  If  two  input  files  are  specified  it  produces  a 
signal  to  noise  ratio  barplot  using  the  first  file  as  the 
numerator  and  the  second  as  the  denominator.  Barplots  are 
produced  on  a  log  scale  so  that  smaller  spikes  will  not  be 
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excluded  front  visibility.  The  rev  value  is  printed  to  the 
left  of  the  plot,  and  in  the  case  of  a  signal  to  noise  plot 
can  be  regarded  as  an  F  ratio  if  the  assumptions  of  section 
3 . 1  are  met. 
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procedure  valeh(var  v:  vactor /start, lgth  v:  integer); 

procedure  to  perform  fast  walah  transform  on  a  vactor  of 
length  n.  Note  that  n  must  be  a  power  of  2.  Also  of 
interaat  ia  that  by  the  symmetry  of  the  FWT  thie  algorithm 
is  lta  own  inverse,  and  can  ba  used  to  generate  valsh 
functions . 


var  temp '.vector? 

half_lgth,  i,  j,  k:  integer; 

begin; 

half_lgth  :»  lgth_v  shr  1; 

j  :*  start; 

for  i  :■  start  to  (start  +  half_lgth  -  1)  do 
begin; 

k  :•  i  +  half  lgth; 
temp [ i ]  v[3)  +  v[j+l); 

temp[k]  :«  v[j]  -  V[j+1); 
if  ((i  -  start)  mod  2)  «  1  then  temp[k) 
j  :  ■=  j  +  2 ; 
end; 

for  i  :■=  start  to  (start  +  lgthv  -  l)  do 
v [ i ]  : «  temp ( i ] ; 

if  (halflgth  >  1)  then 
begin; 

valsh (v, start, half_lgth) ; 
walsh(v, start +half_lgth,half_lgth) ; 
end ; 

end;  (valsh) 


-temp[k] ; 


program  aaq (input, output ) ; 

Program  for  ganaratlng  aaguanclaa  whoaa  croaaterms 
ara  unique.  Generates  permutations,  sorts,  and 
compares  adjacent  elements  for  equality. 

type  vector  •  arrayfl. .210]  of  integer,* 

var  v,p:  vector; 

i,n_param,length,uplim:  integer; 


procedure  permute (var  v:  vector;  n:  integer); 

(generate  the  n  choose  2  interaction  terms  for  n  parameters) 
var  i,j, index;  integer; 

begin; 
index  : *  n ; 

for  i  :■  l  to  (n  *  1)  do 
begin; 

for  j  :«=  (i  +  1)  to  n  do 
begin; 

index  :*  index  +  1; 
v[ index]  :«  v[i]  xor  v[j); 
end; 

end; 

end; (permute) 


procedure  sort (var  v;  vector;  lower, upper;  integer) ; 

(sort  of  a  vector  v  of  length  n  using  recursive  quicksort) 

var  pivot:  integer; 

procedure  partition(var  v:vector ; lower , upper : integer; 

var  pivot : integer) ; 

procedure  swap (var  v: vector; j ,k: integer) ; 

(swaps  the  contents  of  two  vector  locations) 

var  temp:  integer; 

begin; 

temp  : -  v [ j ] ; 
v[j]  v[k] ; 
v[k]  :■  temp; 
end;  (swap)( 


begin; 

pivot  :•  lower; 
repeat 
repeat 

lower  ;■  lower  4  l; 

until  ((v[lower]  >  v[pivot])  or  (lower  >«  upper)); 

while  ((lower  <  upper)  and  (v[pivot]  <  v[upper])) 
do  upper  :«  upper  -  1; 

if  (lower  <  upper)  then  svap(v, lower, upper) ; 

until  (lower  >*  upper) ; 

if  v[ lower]  <  v [pivot]  then 
begin; 

swap(v, lower , pivot) ; 
pivot  :■=  lower; 
end 
else 
begin; 

swap (v, lower-1, pivot) ; 
pivot  :«  lower  -1; 
end; 

end;  (partition) 

begin;  (quicksort) 

part it ion (v, lower, upper, pivot) ; 
if  (pivot  -  lower)  <  (upper  -  pivot)  then 
begin; 

if  ((pivot  -  lower)  >=  2)  then  sort (v, lower, pivot- 
1)  f 

if  ((upper  -  pivot)  >=  2)  then 
sort (v,pivot+l, upper) ; 
end 
else 

begin; 

if  ((upper  -  pivot)  >=  2)  then 
sort (v,pivot+l, upper) ; 

if  ((pivot  -  lower)  >=  2)  then  sort (v, lower, pivot- 

1)  ; 

end; 


end;  (6ort)( 
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function  accapt(var  v:  vector;  n, length:  integer):  boolean 
(accept  or  rajact  a  aat  of  parameter  eeguencies  as  having 
uniqua  alta) 

var  i:  intager; 

result:  boolean; 

begin; 

permute  (v,n)  ; 
sort (v, 1 , length) ; 
result  :■  true; 
i  :«  1; 

repeat 

if  v[i]  *  v[i+l]  then 
result  :«  false; 
i  :■  i  +  1; 

until  (result  «  false)  or  (i  >«  length) ; 
accept  :■«  result; 
end ; { accept } 


procedure  add_lvl(var  v,p:  vector; 
current_lvl ,n_param, tot_lgth:  integer) ; 

{given  unique  elts  at  level  n,  add  elts  at  level  n+1) 
var  i,j, length:  integer; 

begin; 

length  :*=  (sqr (current_lvl)  +  current_lvl)  div  2; 

for  i  :■=  (p[current_lvl-l]+l)  to  (totlgth- (n_param- 
currentlvl) )  do 
begin; 

for  j  :■  1  to  (current_lvl  -  1)  do  v[j]  p [  j ] ; 
v[current_lvl]  :=  i; 
if  accept (v,current_lvl, length)  then 
begin; 

p [ current_l vl ]  : =  i ; 
if  current_lvl  <  n_param  then 

add_lvl (v, p, current_lvl+l, n_param, tot_lgth) 
else 
begin; 

for  j  :■  1  to  current_lvl  do  write (p[ j ] : 4) ; 
writeln; 
end; (write) 
end; (accept  block) 
end; (loop) 
end ; { add  lvl ) { 


write ('enter  A  seguenciee  needed:') 7 
readln (n_param) ; 

length  :■  (eqr (n_param)  +  n_paran)  div  2 
uplin  :■  1; 

while  uplim  <  length  do  uplin  uplin  * 

for  i  :■  1  to  (uplin  -  (njparan  -  1))  do 
begin; 

P[l]  i; 

add_lvl  (v,p,2,n_paran,uplin)  ; 
end; 


mggw ww* ■■jT.'fmufji?!  u  v  .■  gyv? ^,'1  v ? w f 1  ■.■  'i»" y 1  ■.■  -y r*r< 


U*  i'"  t*-  ",*l 
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program  xterm (input, output) ; 

Progran  to  calculata  th*  croaatam  eequency  locations 


type  vector  •  array[1..20]  of  integer; 

var  v:  vector; 

i,n:  integer; 
result;  text; 
fnane:  string [14]; 

procedure  pernute(var  v:  vector;  n;  integer); 

(generate  the  n  choose  2  interaction  terns  for  n  parameters) 
var  i,j, index:  integer; 

begin; 

for  i  ;•  1  to  (n  *  1)  do 

for  j  :*  (i  +  1)  to  n  do 

writeln(result,v[i]  :3, '  xor  ' , v[ j ) : 3 , ' ; ' ,  v[i]  xor 

v[j]:4); 
end; (pernute) 


begin; 

writeln; 

write ('enter  the  output  filename;'); 

read (f name) ; 

assign (result, fname) ; 

rewrite (result) ; 

writeln (result) ; 

write (result, 'enter  the  number  of  sequencies : ' ) ; 
writeln; 

write (' enter  the  number  of  sequencies:'); 

readln(n) ; 

writeln(result,n) ; 

writeln ( 'enter  the  sequencies:'); 

writeln (result, 'enter  the  sequencies: ') ; 

for  i  :*=  1  to  n  do  read (Trm, v[i) )  ; 

for  i  :■  1  to  n  do  write (result, '  ',v[i]); 

writeln (result) ; 

permute (v,n) ; 

close (result) ; 

end. 
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program  spec t rum (input, output ) ; 

{******************************«********»«*****»*************< 
This  program  reads  in  ths  simulation  output  from  a  signal 
run 

and  a  noise  run,  or  from  just  a  signal  run.  It  prompts  the 
user  for  which  type  of  input  he  has,  the  number  of  botches, 
and  size  of  each  batch.  It  than  calculates  the  seguency 
spectrum  if  there  is  only  one  input,  or  calculates  the 
signal 

to  noise  ratio  if  there  are  two.  Plot6  are  on  a  log  scale, 
and  data  values  are  printed  to  the  left  of  the  plot.  Output 
is  written  to  a  file  of  the  users  choosing. 
************************************************************) 

const  size  *  1024; 

half_size  ■  512; 

type  runtype  «  (signal , noise) ; 

vector  ■  array[l. .size]  of  real; 
seguency  «  array [0. .half_size]  of  real; 

var  data:  array [ runtype ]  of  vector; 
run,uplim:  runtype; 
seg:  array [ runtype ]  of  seguency; 
plot:  text; 

infile:  array [runtype]  of  text? 
signal_file,noise_file,plot_file:  string[14] ; 
b_size,n_batches,n_seg, i, j :  integer; { 
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procedure  Initialize; 
begin 

write ('enter  the  name  of  tha  signal  file:'); 

readln(signalfile) ; 

assign (inf ilefsignal ] , signal_f ile)  ; 

reset (infile[signal) ) ; 

writeln; 

writs ( 

'enter  the  name  of  the  noise  file  or  the  word  "nul": 
readln(noise_f ile) ; 
if  (noise_file  «  'nul') 
then 

uplim  :«  signal 
else 
begin 

assign(infile[noise] ,noise_file) ; 
reset (infile[noise) ) ; 
uplim  :«  noise; 
end; 

writeln; 
write ( 

'enter  the  number  of  batches  and  batch  size  of  the 
files: ' )  ; 

readln(n_batches,b_size) ; 
writeln; 

n_seq  :*  b_size  shr  1; 
for  i  :*  0  to  n_seq  do 

for  run  signal  to  uplim  do 
seq[run] [i]  :«  0.0; 

write ('enter  the  name  of  the  output  file:'); 

readln(plot_file) ; 

assign(plot,  plot_file) ; 

rewrite (plot) ; 

writeln; 

write (plot, 'enter  the  name  of  the  signal  file:'); 
writeln (plot, signal_f ile) ; 
writeln (plot) ; 
write (plot, 

'enter  the  name  of  the  noise  file  or  the  word  "nul" 
writeln(plot,noise_file) ; 
writeln (plot) ; 
write (plot, 

'enter  the  number  of  batches  and  batch  size  of  the 
files: ' ) ; 

writeln (plot ,n_batches, '  ' ,b_size) ; 
writeln (plot) ; 
writeln (plot) ; 
end; ( initialize)  { 


procedure  logplot(var  v:  eeguency;  lgth:  integer); 

{A  procedure  to  produce  a  log  acala  barplot  from  vactor  v. 
The  plot  la  acalad  by  tha  largaat  dlffaranca  in  v.) 

var  i,j:  intagar; 

sax, sin, range:  raal; 

begin 

sax  :«  0.0; 
sin  :«  1E37; 
for  i  t»  0  to  lgth  do 
begin 

if  v[i]  >  sax  then  sax  :■  v[i); 

if  (v[i]  <  sin)  and  (v[i)  >  0.0)  then  sin  :■  v[i]; 
end; 

sax  : *  In (sax) ; 
sin  :■  ln(sin/2.8); 
if  sin  >*  sax  then  sin  :«  0; 
range  sax  -  sin; 

for  i  ;*  0  to  lgth  do 
begin  (loop) 

write (plot, v[i] ) ; 
if  (i  sod  5)  -=  0 
then 

vrite(plot,i:5, •+•) 
else 

write  (plot,'  +')»• 
if  v[i]  >  0.0  then 

for  j  :■  l  to  trunc(40* ( (ln(v[i] ) -sin)/range) )  do 
write (plot, •*') ; 
writeln(plot) ; 
end; (loop) 
end; (logplot) 

procedure  crunch (var  v:  vector;  var  seq:  sequency; 

lgth,s_lgth;  integer); 

(procedure  to  calculate  spectrus  estisator  for  vector  v, 
result  placed  in  vector  6eq. ) 

var  i;  integer; 

begin 

seq[0]  seq[0]  +  sqr (v[l]/lgth) ; 
seq[s_lgth]  :«  seq[s_lgth)  +  sqr (v[ lgth )/lgth) ; 

for  i  :«  l  to  6_lgth-l  do 
begin 

seq[i]  :■  seq[i]  +  sqr (v[2*i]/lgth)  + 
sqr (v[  (2*i)+l)/lgth) ; 
end; 

end; (crunch)  { 
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iv' 

A  -.** 


($1  valsh.pas) 
begin 

Initialize; 

for  j  !•  1  to  n_batchee  do 

for  run  :■>  eignal  to  uplin  do 
begin 

for  i  :■  1  to  b_eize  do  read (inf ile [run] ,data [run] [ i] ) ; 
valah (data [run]7l * b_eize) ; 
crunch (data [run] (eeg[run] ,b_size,n_seq) ; 
end; 

if  uplia  «  noiae  then 
for  i  :■  0  to  n_seq  do 

if  (eeq[noise] [ i]  >  0.0)  then 

seq[signal] [i]  :•  seq[signal } [ i ]  /  eeq[noise] [ i] 
else  if ( (seq[noise] [ i]  *  0.0)  and  (seq[ signal ][ i ]  •  0.0)) 

then 

Beq[ signal ][ i]  :«  1.0; 
logplot(seq[ signal] , n_6eq) ; 
close (plot) ; 
end. 


Appendix  2 


Sine*  there  is  no  analytical  method  tor  evaluating  the 
impact  of  a  lag  upon  tha  Walah  spectrum,  va  rasortad  to 
empirical  methods.  It  was  hoped  that  va  would  observe  a 
clear  spike  at  the  original  seguency,  and  at  worst  a  small 
amount  of  noise  elsewhere.  What  we  actually  observed  was 
that  while  the  original  seguency  usually  has  the  largest 
spike,  at  times  it  does  not  have  any  spike  at  all.  We 
present  here  18  spectral  plots  of  different  seguencies  with 
different  lags.  All  plots  were  constructed  from  64 
observations.  The  filename  for  input  i6  descriptive  of  the 
original  seguency  and  lag.  For  example,  LAG602  contains 
Walsh  seguence  number  60  lagged  2  periods.  Recall  that: 
seguency  «  integer [  (seguence  #  +  1)  /  2] 

Notice  that  seguence  numbers  63,  31,  and  15  are 
completely  periodic  with  periods  1,  2  and  4,  respectively.  We 
conseguently  need  not  lag  them  more  than  their  period  lengths 
to  have  completely  identified  their  behavior. 

We  include  the  program  TESTLAG,  which  was  used  to 
generate  the  lagged  Walsh  functions.  The  output  was  plotted 
using  SPECTRUM. 
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program  lagaffact (Input, output) ; 

(MM»*MMMMMM***M*MMM*M*«M**M«***M** 

Genarata  Walsh  functiona  of  a  specified  lag 

A************************************************} 


type  vactor  ■  array [0. .  1023)  of  raal; 

var  v  :  vector; 

seq, size, lag, i  :  Integer; 
outfile  :  text; 
fname  :  string [14]; 

($i  valsh.pas} 

begin; (lagef feet) 

write ('enter  size,  sequence  #,  and  lag:'); 
readln(trm, size, seq, lag) ; 
write('enter  the  output  filename:'); 
readln(trm, fname) ; 
assign (outfile ,  fname); 
rewrite (outfile) ; 

lor  i  :«  0  to  (size  -  1)  do  v[i)  :«=  0; 

'  [6eq]  :■=  1; 

walsh(v, 0,size) ; 

for  i  :■  0  to  (size  -  1)  do 

writeln(outfile,v[ (i-lag)  mod  size):  3); 
close (outfile) ; 
end. (lagef feet) 
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enter  the  naro*  of  the  signal  fils: 
a: lag5l 

sntsr  ths  name  of  the  noise  fils  or  the  word  "nul" : 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


0. 0000000E+00 

Of 

7.8125000E-03 

-(-A****** 

O.OOOOOOOE+OO 

+ 

6 . 6406250E-01 

O.OOOOOOOE+OO 

+ 

3.9062500E-02 

5-|- ******************* 

0.0000000E+00 

+ 

7 . 812500 OE -03 

^-«****** 

O.OOOOOOOE+OO 

+ 

7 . 8125000E-03 

-j-«****** 

O.OOOOOOOE+OO 

10+ 

3 . 9062500E-02 

+******************* 

O.OOOOOOOE+OO 

+ 

3 . 9062500E-02 

+*  ****************** 

O.OOOOOOOE+OO 

+ 

7.8125000E-03 

15-j-******* 

O.OOOOOOOE+OO 

+ 

7.8125000E-03 

+******* 

O.OOOOOOOE+OO 

+ 

3.9062500E-02 

+*  ****************** 

O.OOOOOOOE+OO 

20+ 

3 . 9062500E-02 

+*  ****************** 

O.OOOOOOOE+OO 

+ 

7 . 8125000E-03 

■{-A****** 

O.OOOOOOOE+OO 

+ 

7 . 8125000E-03 

25+«****** 

O.OOOOOOOE+OO 

+ 

3 . 9062500E-02 

-j-*  *****************  * 

O.OOOOOOOE+OO 

+ 

3 . 9062500E-02 

-j-*  ****************** 

O.OOOOOOOE+OO 

30+ 

7 . 8125000E-03 

-j-«****** 

O.OOOOOOOE+OO 

+ 
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enter  the  nan*  of  tha  signal  fila: 
a: lag52 

nter  tha  name  of  tha  noiaa  fila  or  tha  word  "nul": 
nul 

enter  the  number  of  batches  and  batch  aize  of  the  files: 

1  64 


0. OOOOOOOE+OO 
3.1250000E-02 
0. OOOOOOOE+OO 
4 . 0625000E-01 
0. OOOOOOOE+OO 
1.5625000E-01 
0. OOOOOOOE+OO 
3.1250000E-02 
0. OOOOOOOE+OO 
3 . 1250000E-02 
0. OOOOOOOE+OO 
1 . 5625000E-01 
0. OOOOOOOE+OO 
1.5625000E-01 
0. OOOOOOOE+OO 
3.1250000E-02 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 


0+ 

+**  ********* 

+ 

+ 

+ 

+*********** 

+ 

+*********** 

1°+ 

+ 

+ 

15+*********** 

+ 

+ 

+ 

+ 

20+ 

+ 

+ 

+ 

+ 

25+ 

+ 

+ 

+ 

+ 

30+ 

+ 

+ 


59 


enter  the  nan*  of  tha  signal  fila: 
a : lag53 

antar  tha  name  of  tha  noise  fila  or  tha  word  HnulM: 
nul 

enter  tha  number  of  batches  and  batch  size  of  the  files: 
1  64 


0. OOOOOOOE+OO 

o+ 

7.0312500E-02 

+************************** 

0. OOOOOOOE+OO 

+ 

2.2656250E-01 

0. OOOOOOOE+OO 

+ 

3.5156250E-01 

5+*  ************************  * 

0. OOOOOOOE+OO 

+ 

7 . 0312500E-02 

+************************** 

0. OOOOOOOE+OO 

+ 

7.8125000E-03 

X******** 

0. OOOOOOOE+OO 

10+ 

3.9062500E-02 

+*  ******************** 

0. OOOOOOOE+OO 

+ 

3.9062500E-02 

-|-********************* 

0. OOOOOOOE+OO 

+ 

7.8125000E-03 

15-j-******** 

0. OOOOOOOE+OO 

+ 

7.8125000E-03 

+******** 

0. OOOOOOOE+OO 

+ 

3 .9062500E-02 

0.  OOOOOOOE+OO 

20+ 

3.9062500E-02 

^-********************* 

0.  OOOOOOOE+OO 

+ 

7 . 8125000E-03 

-|-******** 

0. OOOOOOOE+OO 

+ 

7.8125000E-03 

25+******** 

0. OOOOOOOE+OO 

+ 

3 . 9062500E-02 

0. OOOOOOOE+OO 

+ 

3.9062500E-02 

0. OOOOOOOE+OO 

30+ 

7.8125000E-03 

-j-******** 

0. OOOOOOOE+OO 

+ 

enter  the  name  of  the  signal  file: 
a: lag54 

enter  the  name  of  the  noise  file  or  the  word  MnulM: 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


0. OOOOOOOE+OO 
1 . 2500000E-01 
0. 0000000E+00 
1.2500000E-01 
0. OOOOOOOE+OO 
6. 2500000E-01 
0. OOOOOOOE+OO 
1.2500000E-01 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0. OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0.  OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 
0. OOOOOOOE+OO 


+*************** 

+ 

+ 

+ 

+ 

10+ 

+ 

t 

+ 

+ 

15+ 

+ 

+ 

t 

+ 

20+ 

+ 

t 

+ 

+ 

25+ 

t 

+ 

+ 

+ 

30+ 

+ 

+ 
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enter  the  name  of  the  signal  file: 
a: lag55 

enter  the  name  of  the  nolee  file  or  the  word  "nul" : 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


O.OOOOOOOE-t-OO  Of 

1.9531250E-01  +* ********  ************************** 

O.OOOOOOOE+OO  + 

1 . 0156250E-01  f* **************** ************ 

O.OOOOOOOE+OO  + 

3.5156250E-01  5f** ******************** ****************** 

O.OOOOOOOE+OO  f 

7 . 03 12  500E-02  +** *********** ************* 

O.OOOOOOOE+OO  + 

7 . 8 125000E-03  f ******** 

O.OOOOOOOE+OO  10f 

3 . 9062500E-02  t***** **************** 

O.OOOOOOOE+OO  f 

3 . 9062500E-02  f******* ************** 

O.OOOOOOOE+OO  f 

V.8125000E-03  15f+******* 

O.OOOOOOOE+OO  f 

7 . 8 125000E-03  f******** 

O.OOOOOOOE+OO  f 

3.9062500E-02  f ********************* 

O.OOOOOOOE+OO  20f 

3 . 9062500E-02  f********** *********** 

O.OOOOOOOE+OO  -j- 

7 . 8125000E-03  f******** 

O.OOOOOOOE+OO  f 

7 . 8125000E-03  25f******** 

O.OOOOOOOE+OO  -j- 

3 . 9062500E-02  -f* ******************** 

O.OOOOOOOE+OO  f 

3 . 9062500E-02  t********************* 

O.OOOOOOOE+OO  30+ 

7 . 8125000E-03  +******** 

O.OOOOOOOE+OO  f 


o 


enter  the  nan*  of  the  signal  fils: 
a : lag60l 

enter  the  name  of  the  noiae  file  or  the  word  HnulM: 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


0. OOOOOOOE+OO 

O.OOOOOOOE+OO 

3.1250000E-02 

0.0000000E+00 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

3 . 12  50000E-02 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

3 . 12  50000E-02 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

3 . 1250000E-02 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

3 . 1250000E-02 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

3.1250000E-02 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

3 . 1250000E-02 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

7.8125000E-01 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 


0+ 

4 

■j-AAAAAAAAA 

+ 

+ 

5+ 

-!********* 

+ 

+ 

4 

10+*  ******** 

+ 

t 

4 

+AAAAAAAAA 

15+ 

+ 

4 

-(-A******** 

+ 

20+ 

4 

+AAAA AAAAA 

+ 

4 

2  5+ 

-j-AAAAAAAAA 

+ 

+ 

+ 

+ 

+ 


.  *•  '■  ■.***',« 
M..  uft*.  .*A*. L-'m  V  ... 
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enter  the  name  of  the  signal  fils: 
a : lag602 

enter  the  name  of  the  noise  file  or  the  word  "nul": 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


0.0000000E+00 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

1 . 2500000E-01 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

1 . 2500000E-01 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

1 . 2500000E-01 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

6.2500000E-01 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 


0+ 

+ 

+ 

+ 

+ 

5+ 

+ 

+ 

+ 

+ 

10+ 

+ 

+ 

+ 

+ 

15+ 

+ 

+ 

************** 

+ 

20+ 

+ 

-}-*************** 

+ 

+ 

2  5+ 

+ 

+ 

+ 

+ 

+ 
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enter  the  name  of  tha  signal  fils: 

a:lag603 

sntsr  ths  nans  of  ths  noise  file  or  the  word  Hnul*': 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 

1  64 


O.OOOOOOOE+OO  0+ 

O.OOOOOOOE+OO  + 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  5+ 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  + 

0.0000000E+00  + 

O.OOOOOOOE+OO  + 

3.1250000E-02  10-j-********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  15+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  20+ 

O.OOOOOOOE+OO  + 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  25+ 

2 . 8125000E-01  +********************************* 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

5.3125000E-01  30+**** ************************************ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 


enter  the  name  of  the  signal  fils: 
a : lag604 


sntsr  ths  name  of  ths  nolss  fils  or  ths  word  "nul": 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


O.OOOOOOOE+OO  0+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  5+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  10+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  15+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  20+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  25+ 

5. 0000000E-01  +**************************************** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

5 . OOOOOOOE-Ol  30+* *************************************** 
O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 
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enter  the  nan*  of  tha  signal  fila: 
a: lag605 

antar  tha  name  of  tha  noise  fila  or  tha  word  "nul" : 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


O.OOOOOOOE+OO  0+ 

O.OOOOOOOE+OO  + 

3.1250000E-02  +********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  5+ 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

3 . 1250000E-02  10+********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  15+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  20+ 

O.OOOOOOOE+OO  + 

3 . 1250000E-02  +********** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  25+ 

2 . 8125000E-01  +********************************* 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

5 . 3125000E-01  30+**************************************** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 


s  vs 
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enter  the  nan*  of  the  signal  file: 
a : lag608 


-  A 

fj 


enter  the  name  of  the  noiee  file  or  the  word  "nul": 
nul 


enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


0. OOOOOOOE+OO 

O.OOOOOOOE+OO 

0. OOOOOOOE+OO 

O.OOOOOOOE+OO 

0. 0000000E+00 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

1. OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 


0+ 

+ 

+ 

+ 

+ 

5+ 

+ 

+ 

+ 

+ 

10+ 

+ 

+ 

+ 

t 

15+ 

+ 

+ 

+ 

+ 

20+ 

+ 

+ 

+ 

+ 

25+ 

+ 

+ 

+ 

+ 

+ 

+ 


enter  the  nan*  of  tha  signal  flla 
a : lag6016 


enter  tha  name  of  tha  noise  fila  or  tha  word  Mnulw: 
nul 

antar  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

1 . OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 


0+ 

+ 

+ 

+ 

+ 

5+ 

t 

+ 

t 

+ 

10+ 

+ 

+ 

+ 

+ 

15+ 

+ 

+ 

+ 

+ 

20+ 

+ 

+ 

+ 

+ 

25+ 

+ 

+ 

+ 

+ 

+ 

+ 
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enter  the  nan*  of  tha  signal  file: 
a:lag631 

enter  tha  nana  of  tha  nolaa  fila  or  the  word  "nul": 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


0. 0000000E+00  0+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  5+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  10+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  15+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  20+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  25+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  30+ 

O.OOOOOOOE+OO  + 

1. OOOOOOOE+OO  +**************************************** 


->  ^  V-  .  V  " 
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enter  the  name  of  tha  signal  fila: 
a:lag311 

antar  tha  name  of  tha  noiaa  fila  or  tha  word  "nul": 
nul 

antar  tha  number  of  batches  and  batch  size  of  the  files: 
1  64 


0.0000000E+00  Of 
O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  5+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  10+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  15+ 

1. 0000000E+00  +**************************************** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  20+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  25+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  30+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 


enter  the  name  of  the  signal  fils: 
a : lag312 

enter  the  name  of  the  noise  file  or  the  word  MnulM : 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


O.OOOOOOOE+OO  0+ 

0. OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  5+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

0. OOOOOOOE+OO  + 

0.0000000E+00  10+ 

0.0000000E+00  + 

0.0000000E+00  + 

0.0000000E+00  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  15+ 

1 .  OOOOOOOE+OO  +**************************************** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  20+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  25+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  30+ 

0. OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 
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enter  the  nan*  of  th*  signal  fils: 
a : lagl51 

•nt*r  th*  nan*  of  th*  nois*  fil*  or  th*  word  Mnul" : 
nul 

•nt*r  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


O.OOOOOOOE+OO  0+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

0.0000000E+00  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  5+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

5.0000000E-01  |**************************************** 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  10+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  15+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  20+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

5 . OOOOOOOE-Ol  +**************************************** 

O.OOOOOOOE+OO  25+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  30+ 

O.OOOOOOOE+OO  + 

O.OOOOOOOE+OO  + 
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enter  the  name  of  the  signal  file: 
a : lagl52 

enter  the  name  of  the  noise  file  or  the  word  ”nul": 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

1 . OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 


Of 

+ 

+ 

+ 

i 

5+ 

+ 

+ 

+ 

10+ 

+ 

+ 

+ 

t 

15+ 

+ 

t 

+ 

+ 

20+ 

t 

+ 

+ 

+ 

25+ 

+ 

+ 

+ 

+ 

30+ 

+ 

t 
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enter  the  name  of  the  signal  file: 
a: lagl53 

enter  the  name  of  the  noiee  file  or  the  vord  "nul” : 
nul 

enter  the  number  of  batches  and  batch  size  of  the  files: 
1  64 


O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

5 . OOOOOOOE-Ol 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

5. OOOOOOOE-Ol 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 


0+ 

+ 

+ 

+ 

+ 

5+ 

t 

t 

+ 

10+ 

+ 

+ 

+ 

+ 

15+ 

+ 

+ 

+ 

+ 

20+ 

+ 

+ 

+ 

25+ 

+ 

+ 

+ 

+ 

30+ 

+ 
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present  a  promising  method  for  identifying  significant 
parameters.  However,  the  method  exhibits  undesirable  behavior 
when  time  lag6  are  present  in  the  model. 


